#*********************************************************#
#                                                         #
#                                                         #
#                 REPLICATION FILES FOR:                  #
#                                                         #
# Batons and Ballots: The Effectiveness of State Violence #
#        in Fighting Against Catalan Separatism           #                                       
#                                                         #         
#                                                         #
#                                                         #
#*********************************************************#
#
#***************#
# Load Packages #
#***************#

require(MatchIt)
require(reshape2)
require(ggplot2)
require(Rmisc)
require(Hmisc)

#############
##Load data##
#############
setwd('file_path')
X <- read.csv2("Replication_Data_RAP.csv")

###############
####TABLE 1####
###############

#Table 1: Column 1
summary(m1 <- lm(Independence2017 ~ Violence + Independence2015, data = X))

#Table 1: Column 2
summary(m2 <- lm(Independence2017 ~ Violence + Independence2015 + Turnout2015 + log(Census2015) + log(cap_dist+0.1) + as.factor(Province), data = X))

#Matching algorithm

all.out_cem <- matchit(Violence ~ as.numeric(as.character(Independence2015)) + as.numeric(as.character(Turnout2015)) + as.numeric(as.character(Census2015)) + as.numeric(as.character(cap_dist)) + as.factor(Province), data = X[!X$Town == 'Girona' & !X$Town == 'Barcelona' & !X$Town == 'Lleida' & !X$Town == 'Tarragona',], method = 'cem')

matchdata_cem <- match.data(all.out_cem)

#Table 1: Column 3

summary(mat1_cem <- lm(Independence2017 ~ Violence + Independence2015, data = matchdata_cem, weights = weights))

#Table 1: Column 4

summary(mat2_cem <- lm(Independence2017 ~ Violence + Independence2015 
                       + Turnout2015 + log(Census2015) 
                       + log(cap_dist+0.1) + as.factor(Province), 
                       data = matchdata_cem, 
                       weights = weights
                       ))


##Stack Data for DID

d0 <- reshape(X, direction="long", varying=list(names(X)[7:8]), v.names=c("Independence"), idvar=c("Town"), timevar="Year", times=c(2015, 2017))

paneldata0 <- as.data.frame(cbind(d0$Town, d0$Year, d0$Independence, ifelse(d0$Year == 2015, d0$Turnout2015, d0$Turnout2017), ifelse(d0$Year == 2015, d0$Census2015, d0$Census2017), d0$Violence, d0$Violence_H, d0$Violence_L))

colnames(paneldata0) <- c('Town', 'Year', 'Independence', 'Turnout', 'Census', 'Violence', 'Violence_H', 'Violence_L')

paneldata0$Independence <- as.numeric(as.character(paneldata0$Independence))
paneldata0$Census <- as.numeric(as.character(paneldata0$Census))
paneldata0$Turnout <- as.numeric(as.character(paneldata0$Turnout))

#Table 1: Column 5
summary(DID0 <- lm(Independence ~ as.factor(Year)*Violence, data = paneldata0))

#Table 1: Column 6
summary(DID1 <- lm(Independence ~ as.factor(Year)*Violence + log(Census) + Turnout, data = paneldata0))

##Matched Stack Data for DID-M

d <- reshape(matchdata_cem, direction="long", varying=list(names(matchdata_cem)[7:8]), v.names=c("Independence"), idvar=c("Town"), timevar="Year", times=c(2015, 2017))

paneldata_m0 <- as.data.frame(cbind(d$Town, d$Year, d$Independence, ifelse(d$Year == 2015, d$Turnout2015, d$Turnout2017), ifelse(d$Year == 2015, d$Census2015, d$Census2017), d$Violence, d$Violence_H, d$Violence_L, d$weights))

colnames(paneldata_m0) <- c('Town', 'Year', 'Independence', 'Turnout', 'Census', 'Violence', 'Violence_H', 'Violence_L', 'Weights')

paneldata_m0$Independence <- as.numeric(as.character(paneldata_m0$Independence))
paneldata_m0$Census <- as.numeric(as.character(paneldata_m0$Census))
paneldata_m0$Turnout <- as.numeric(as.character(paneldata_m0$Turnout))

#Table 1: Column 7
summary(DIDM0 <- lm(Independence ~ as.factor(Year)*Violence, data = paneldata_m0, weights = Weights))

#Table 1: Column 8
summary(DIDM1 <- lm(Independence ~ as.factor(Year)*Violence + log(Census) + Turnout, data = paneldata_m0, weights = Weights))


#Figure 2a
tgc0 <- summarySE(paneldata0, measurevar="Independence", groupvars=c("Violence","Year"))
tgc0

tgc0$Violence <- as.factor(tgc0$Violence)
tgc0$Year <- as.factor(as.character(tgc0$Year))

pd <- position_dodge(0.1)

full_ind <- ggplot(tgc0, aes(x=Year, y=Independence,
                             colour=Violence, group=Violence)) + 
  geom_errorbar(aes(ymin=Independence-se, ymax=Independence+se), width=.1, position=pd) +
  geom_line(position=pd, size = 1) +
  geom_point(position=pd, size=3, shape=19) +
  xlab("Year") +
  ylab("Support for Pro-Independence Parties (Percentage)") +
  theme_bw() +
  scale_color_manual(labels = c("No", "Yes"), values = c("blue", "red")) +
  labs(color='State Violence') + #Set legend title
  theme(legend.justification=c(1,0),
        legend.position=c(0.995,0.01),
        axis.text=element_text(size=18),
        axis.title=element_text(size=22,face="bold"),
        legend.text=element_text(size=14))

pdf("Plot_Full_Independence.pdf", height = 6, width = 7)
full_ind
dev.off()

#Figure 2b
tgc <- summarySE(paneldata_m0, measurevar="Independence", groupvars=c("Violence","Year"))[,1:6]
tgc

a <- paneldata_m0[which(paneldata_m0$Violence == 0 & paneldata_m0$Year == 2015),]
b <- paneldata_m0[which(paneldata_m0$Violence == 0 & paneldata_m0$Year == 2017),]
c <- paneldata_m0[which(paneldata_m0$Violence == 1 & paneldata_m0$Year == 2015),]
d <- paneldata_m0[which(paneldata_m0$Violence == 1 & paneldata_m0$Year == 2017),]

tgc[1,3] <- length(a$Independence)
tgc[1,4] <- wtd.mean(a$Independence, a$Weights)
tgc[1,5] <- sqrt(wtd.var(a$Independence, a$Weights))
tgc[1,6] <- tgc[1,5] / sqrt(tgc[1,3])

tgc[2,3] <- length(b$Independence)
tgc[2,4] <- wtd.mean(b$Independence, b$Weights)
tgc[2,5] <- sqrt(wtd.var(b$Independence, b$Weights))
tgc[2,6] <- tgc[2,5] / sqrt(tgc[2,3])

tgc[3,3] <- length(c$Independence)
tgc[3,4] <- wtd.mean(c$Independence, c$Weights)
tgc[3,5] <- sqrt(wtd.var(c$Independence, c$Weights))
tgc[3,6] <- tgc[3,5] / sqrt(tgc[3,3])

tgc[4,3] <- length(d$Independence)
tgc[4,4] <- wtd.mean(d$Independence, d$Weights)
tgc[4,5] <- sqrt(wtd.var(d$Independence, d$Weights))
tgc[4,6] <- tgc[4,5] / sqrt(tgc[4,3])

tgc$Violence <- as.factor(tgc$Violence)
tgc$Year <- as.factor(as.character(tgc$Year))
pd <- position_dodge(0.1)

match_ind <- ggplot(tgc, aes(x=Year, y=Independence,
                             colour=Violence, group=Violence)) + 
  geom_errorbar(aes(ymin=Independence-se, ymax=Independence+se), width=.1, position=pd) +
  geom_line(position=pd, size = 1) +
  geom_point(position=pd, size=3, shape=19) +
  xlab("Year") +
  ylab("Support for Pro-Independence Parties (Percentage)") +
  theme_bw() +
  scale_color_manual(labels = c("No", "Yes"), values = c("blue", "red")) +
  labs(color='State Violence') + #Set legend title
  theme(legend.justification=c(1,0),
        legend.position=c(0.995,0.01),
        axis.text=element_text(size=18),
        axis.title=element_text(size=22,face="bold"),
        legend.text=element_text(size=14)) # Position legend in bottom right

pdf("Plot_Match_Independence.pdf", height = 6, width = 7)
match_ind
dev.off()

###############
####TABLE 2####
###############

#Table 2: Column 1
summary(m3 <- lm(Turnout2017 ~ Violence + Turnout2015, data = X))

#Table 2: Column 2
summary(m4 <- lm(Turnout2017 ~ Violence + Turnout2015 + Independence2015 + log(Census2015) + log(cap_dist+0.1) + as.factor(Province), data = X))

#Table 2: Column 3
summary(mat3 <- lm(Turnout2017 ~ Violence + Turnout2015, data = matchdata_cem, weights = weights))

#Table 2: Column 4
summary(mat4 <- lm(Turnout2017 ~ Violence + Turnout2015 + Independence2015 + log(Census2015) + log(cap_dist+0.1) + as.factor(Province), data = matchdata_cem, weights = weights))

#Table 2: Column 5
summary(DIDM0 <- lm(Turnout ~ as.factor(Year)*Violence, data = paneldata0))

#Table 2: Column 6
summary(DIDM1 <- lm(Turnout ~ as.factor(Year)*Violence + log(Census) + Independence, data = paneldata0))

#Table 2: Column 7
summary(DIDM0 <- lm(Turnout ~ as.factor(Year)*Violence, data = paneldata_m0, weights = Weights))

#Table 2: Column 8
summary(DIDM1 <- lm(Turnout ~ as.factor(Year)*Violence + log(Census) + Independence, data = paneldata_m0, weights = Weights))

#Figure 2c

tgc1 <- summarySE(paneldata0, measurevar="Turnout", groupvars=c("Violence","Year"))
tgc1

tgc1$Violence <- as.factor(tgc1$Violence)
tgc1$Year <- as.factor(as.character(tgc1$Year))

pd <- position_dodge(0.1)

full_turn <- ggplot(tgc1, aes(x=Year, y=Turnout,
                             colour=Violence, group=Violence)) + 
  geom_errorbar(aes(ymin=Turnout-se, ymax=Turnout+se), width=.1, position=pd) +
  geom_line(position=pd, size = 1) +
  geom_point(position=pd, size=3, shape=19) +
  xlab("Year") +
  ylab("Turnout") +
  theme_bw() +
  scale_color_manual(labels = c("No", "Yes"), values = c("blue", "red")) +
  labs(color='State Violence') + #Set legend title
  theme(legend.justification=c(1,0),
        legend.position=c(0.995,0.01),
        axis.text=element_text(size=18),
        axis.title=element_text(size=22,face="bold"),
        legend.text=element_text(size=14))

pdf("Plot_Full_Turnout.pdf", height = 6, width = 7)
full_turn
dev.off()

#Figure 2d

tgc2 <- summarySE(paneldata_m0, measurevar="Turnout", groupvars=c("Violence","Year"))[,1:6]
tgc2

require(Hmisc)
tgc2[1,3] <- length(a$Turnout)
tgc2[1,4] <- wtd.mean(a$Turnout, a$Weights)
tgc2[1,5] <- sqrt(wtd.var(a$Turnout, a$Weights))
tgc2[1,6] <- tgc2[1,5] / sqrt(tgc2[1,3])

tgc2[2,3] <- length(b$Turnout)
tgc2[2,4] <- wtd.mean(b$Turnout, b$Weights)
tgc2[2,5] <- sqrt(wtd.var(b$Turnout, b$Weights))
tgc2[2,6] <- tgc2[2,5] / sqrt(tgc2[2,3])

tgc2[3,3] <- length(c$Turnout)
tgc2[3,4] <- wtd.mean(c$Turnout, c$Weights)
tgc2[3,5] <- sqrt(wtd.var(c$Turnout, c$Weights))
tgc2[3,6] <- tgc2[3,5] / sqrt(tgc2[3,3])

tgc2[4,3] <- length(d$Turnout)
tgc2[4,4] <- wtd.mean(d$Turnout, d$Weights)
tgc2[4,5] <- sqrt(wtd.var(d$Turnout, d$Weights))
tgc2[4,6] <- tgc2[4,5] / sqrt(tgc2[4,3])

tgc2$Violence <- as.factor(tgc2$Violence)
tgc2$Year <- as.factor(as.character(tgc2$Year))

pd <- position_dodge(0.1)

match_turn <- ggplot(tgc2, aes(x=Year, y=Turnout,
                              colour=Violence, group=Violence)) + 
  geom_errorbar(aes(ymin=Turnout-se, ymax=Turnout+se), width=.1, position=pd) +
  geom_line(position=pd, size = 1) +
  geom_point(position=pd, size=3, shape=19) +
  xlab("Year") +
  ylab("Turnout") +
  theme_bw() +
  scale_color_manual(labels = c("No", "Yes"), values = c("blue", "red")) +
  labs(color='State Violence') + #Set legend title
  theme(legend.justification=c(1,0),
        legend.position=c(0.995,0.01),
        axis.text=element_text(size=18),
        axis.title=element_text(size=22,face="bold"),
        legend.text=element_text(size=14))

pdf("Plot_Match_Turnout.pdf", height = 6, width = 7)
match_turn
dev.off()
